	   	    
     
      Program Moisture
      
	parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
	implicit double precision (a-h,o-z)
	common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	character chartin*1,infile*20,outfile*20
	common/comc/chartin,infile,outfile
	common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping
	npppsi = 4
c	in outp, for each material layer print out values at npppsi+1 points
      call initl
	call setup
	call outp(npppsi)
	do 100 index = 1, ntsteps
	   kdt = index
	   t = t + dt
	   call advance
	   iscr = kdt / levt
	   iscr = iscr * levt
	   tout = t/toutfac
c	               if(kdt  .eq.  (levt/20))  call outp(npppsi)
c 	   above call to output to monitor ``short time`` behavior of numerical
c	   solution
	   if(iscr  .eq.  kdt)  call outp(npppsi)
 100   continue
	close(66)
	close(67)
      stop 'Fin du Programme de diffuse'											   
	 end
c=======================================================
c=======================================================
	subroutine advance
c	   advance the solution one dt time step
c	   have matrix al, ad, ar of tridiagonal linear system
c	   form right hand side of current (nth time step) solution
c	   present value of t is at ``new`` time level (n+1 time level) where
c	   are about to calculate solution values
     
      parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
	  implicit double precision (a-h,o-z)
 	  common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	  character chartin*1,infile*20,outfile*20
	  common/comc/chartin,infile,outfile
	  common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping

	onema = one - alpha
	three = 3.d0
	six = 6.d0
	k = 1
	x = xpta(k)
	told = t - dt
	ubdyn = dirl(told)
	ubdynpl = dirl(t)
	if(initu  .eq.  1) ubdyn = truesol(zero, told)
	if(initu  .eq.  1) ubdynpl = truesol(zero, t)
	rtside(k) = aln(k)*ubdyn + adn(k)*ua(k) + arn(k)*ua(k+1)
	rtside(k) = rtside(k) - al(k)*ubdynpl
	if(initu  .eq.  1)  then
	   xl = xpta(k-1)
	   xr = xpta(k+1)
	   hk = dxa(k)
	   hkp = dxa(k+1)
	   hk6 = hk/six
	   hk3 = hk/three
	   hkp6 = hkp/six
	   hkp3 = hkp/three
c	   level n+1 terms
	   sl = source(xl,t)
	   sd = source(x,t)
	   sr = source(xr,t)
	   rtside(k) = rtside(k)+alpha*(sl*hk6 + sd*(hk3+hkp3) + sr*hkp6)
c	   level n terms
	   sl = source(xl,told)
	   sd = source(x,told)
	   sr = source(xr,told)
	   rtside(k) = rtside(k)+onema*(sl*hk6 + sd*(hk3+hkp3) + sr*hkp6)
	endif
c	
	nunknm1 = nunkn - 1
	do  k = 2, nunknm1
	   x = xpta(k)
	   rtside(k) = aln(k)*ua(k-1) + adn(k)*ua(k) + arn(k)*ua(k+1)
	   if(initu  .eq.  1)  then
	      xl = xpta(k-1)
	      xr = xpta(k+1)
	      hk = dxa(k)
	      hkp = dxa(k+1)
	      hk6 = hk/six
	      hk3 = hk/three
	      hkp6 = hkp/six
	      hkp3 = hkp/three
c	      level n+1 terms
	      sl = source(xl,t)
	      sd = source(x,t)
	      sr = source(xr,t)
	 rtside(k) =rtside(k)+alpha*(sl*hk6 + sd*(hk3+hkp3)+ sr*hkp6)
c	      level n terms
	      sl = source(xl,told)
	      sd = source(x,told)
	      sr = source(xr,told)
	 rtside(k)=rtside(k)+onema*(sl*hk6 + sd*(hk3+hkp3)+ sr*hkp6)
	   endif
      enddo
c	
	k = nunkn
      x = xpta(k)
	ubdyn = dirl(told)
	ubdynpr = dirl(t)
	if(initu  .eq.  1) ubdyn = truesol(xrt, told)
	if(initu  .eq.  1) ubdynpr = truesol(xrt, t)
	rtside(k) = aln(k)*ua(k-1) + adn(k)*ua(k) + arn(k)*ubdyn
	rtside(k) = rtside(k) - ar(k)*ubdynpr
	   if(initu  .eq.  1)  then
	      xl = xpta(k-1)
	      xr = xpta(k+1)
	      hk = dxa(k)
	      hkp = dxa(k+1)
	      hk6 = hk/six
	      hk3 = hk/three
	      hkp6 = hkp/six
	      hkp3 = hkp/three
c	      level n+1 terms
	      sl = source(xl,t)
	      sd = source(x,t)
	      sr = source(xr,t)
	 rtside(k)= rtside(k)+alpha*(sl*hk6 + sd*(hk3+hkp3) + sr*hkp6)
c	      level n terms
	      sl = source(xl,told)
	      sd = source(x,told)
	      sr = source(xr,told)
	rtside(k)=rtside(k)+onema*(sl*hk6 + sd*(hk3+hkp3) + sr*hkp6)
	   endif
c	
	call trisol(al,ad,ar,nunkn,rtside,ua(1))
	ua(0) = ubdynpl
	if(isymrt  .eq.  0) ua(nunkn+1) = ubdynpr
	 return
	 end
c=======================================================
c=======================================================
	 function dirl(time)
	parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
	  implicit double precision (a-h,o-z)
	  common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	 character chartin*1,infile*20,outfile*20
	common/comc/chartin,infile,outfile
	common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping
	 dirl = outrhl / 80.d0
	 return
c	
	 entry dirr(time)
	 dirr = outrhr / 80.d0
	 return
	 end
c=======================================================
c=======================================================
	 function truesol(xx, time)
	 implicit double precision(a-h,o-z)
c	 truesol = 8.d0 * time * time + 7.d0 * xx * time +.0 * xx * xx
	 pid2 = 1.57079632679d0
	 if(xx  .lt.  pid2)  then
	     truesol = sin(xx) * exp(-50.d0*time)
	 else 
	     truesol = sin(5.0d0*xx) * exp(-50.d0*time)
	           end if
	           return
c	
	entry source(xx,time)
c	source = 16.d0 * time + 7.d0 * xx - 30.d0
c	source = source * .81d0
	source = 0.d0
	return
	end
c===========================================================
c===========================================================
cc    FORTRAN SOURCE FILE : initl.for
	subroutine initl
c	read in data from input file (physical constants and specifications)
c      For the finite difference scheme)
	 parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
	implicit double precision (a-h,o-z)
	common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	 character chartin*1,infile*20,outfile*20
	common/comc/chartin,infile,outfile
	common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping
	 character comment*72
	 character cscr*72
	 character outcol*72
	 zero = 0.d0
	 one = 1.d0
	 two = 2.d0
	 write(*,*)'type in file name of the input data'
	 write(*,*)'(up to 20 characters  do NOT encloe in apostrophes)'
	 write(*,*)'***LIMIT length of name to machine limit'
	 write(*,*)'***do NOT use any leading blanks'
	 read(*,500) infile
500         format(a20)
	 write(*,*)'type in file name for the output'
	 write(*,*)'(up to 20 characters  do NOT encloe in apostrophes)'
	 write(*,*)'***LIMIT length of name to machine limit - 3'
	 write(*,*)'***do NOT use any leading blanks'
	 write(*,*)'note columnformat output file is also created'
	 write(*,*)'its file name = <output file name>col'
	 read(*,500) outfile
	 open(55, file=infile)
	 open(66, file=outfile)
	 write(66,*)'  output file name is  ' ,outfile
	 write(66,*)'  input file name is  ' ,infile
c	get outcol name
	   if(outfile(1:1) .eq.'  ')  then
	write(*,*)'WARNINGfirst char of outfile is blank,outcol=col'
	      outcol = 'col'
	      goto  11
	   endif
	do index = 2,20
	   if(outfile(index : index)  .eq.  '  ')  then
	      outcol = outfile(1 : index-1)//'col'
	      goto  11
	   endif
	enddo
	outcol = outfile//'col'
11    continue
	open(67,file=outcol)
	write(67,*)'  column output file name is  ' ,outcol
	write(67,*)'  input file name is  ' ,infile
	read(55,*) nlayers
	write(66,*)'  number of material layers is  ' ,nlayers
c	check dimensions
	if(nlayers.gt.nlmax) stop 'nlayers gt nlmax   increase nlmax'
	read(55,*) rhsol
	write(66,*)' rhsol = relative humidity in % at which the' ,nlayers
	write(66,*)'  equilibrium solubility (wt  %) of H2O in all the '
	write(66,661)'material layers was measured;rhsol  (%) = ' ,  rhsol
661   format(a, f6.2)
	rhfac = 80.d0/rhsol
c	use rhfac to convert everything to case where rhsol = 80% as is
c	case for original coe
	         do 10  k=1,nlayers
	            write(66,*)  '             layer number  ',k
	 read(55,*) thick
	 read(55,*) diffcf
	 read(55,*) rho
	 read(55,*) solwp
	 thicka(k) = thick
	 diffcfa(k) =  diffcf
	 rhoa(k) = rho
	 solwpa(k) = solwp
	              write(66,*)  thicka(k),  '  thicka(cm)'
	              write(66,*)  diffcfa(k),  '  diffcfa(cm*cm/sec)'
	              write(66,*)  rhoa(k),  '  rhoa(g/cc)'
	
     
      write(66,*)  solwpa(k), ' solwpa(wt% @rhsol%RH gH2O/100g dry m) '
c	convert to case where rhsol is 80%
	 solwpa(k) = solwp * rhfac
10           continue
	         isymrt = 0
	write(66,*) isymrt,'  if  0  (1) right end is exposed (sealed)'
	         kleftbc = 0
c	Dirichlet boundary data at left endpoint (solution value specified)
	         krtbc = isymrt
c	krbtc = 0 for Dirichlet data at right endpoint,
c	krbtc = 1 for symmetry condition (0 Neumann data)  at right endpoint
	read(55,*) outrhl, outrhr
	write(66,*) outrhl, outrhr,'outside relative humidity(%)'
	write(66,*) 'atthe left,right bound(endofthe material)'
	          npscr = 0
	          do 20  k=1,nlayers
	          read(55,*) nintlr
	  nintlra(k) = nintlr
	write(66,*)  nintlra(k),'# finite diff subintervals in layer',k
	  npscr = npscr + nintlr
20            continue
	  npscr = npscr + 1
c	check dimensions
	if(npscr  .gt.  nptsmax) stop'npts gt ntsmax increase nptsmax'
	          read(55,*) chartin
	          cscr =  '  '//chartin
	write(66,*) cscr,' char var of length  1  ---  input time units'
	write(66,*)' time units are: s = seconds,  i = minutes,  '
	write(66,*)'time units are: h = hours,  d = days,  w = weeks,  '
	write(66,*)'time units are:m= months  (30.43667  days/month)'
	write(66,*)'time units are: y = years  (365.24  days/year)'
	          read(55,*) tf
	write(66,*)tf,'final time (in given units)  calc carried out to'
	          read(55,*) ntsteps
	write(66,*) ntsteps,'#  of finite diff time  steps to get to tf'
	          read(55,*) levt
	write(66,*)levt,'print approx soln values every levt time steps'
	          read(55,*) initu
	write(66,*) initu,'if 0 (1)  start with 0 (known test)  concen. '
	          read(55,*) alpha
	write(66,*)alpha, 'time ave parameter,1. = implicit, .5  =  C.N. '
	          read(55,*) lumping
	write(66,*)lumping, '  if 1 lump mass matrix, else do not lump'
c	any remaining lines (up to 5) in input file are comments to be
c	printed in the output file
	          do 30 k=1,5
	           read(55,501,end=31) comment
501           format(a72)
	           write(66,601) comment
601           format('  ',a72)
30             continue
31             continue
	 close(55)
	 return
	 end
c=======================================================
c=======================================================
	 subroutine tridec(al,ad,ar,nrows)
	 implicit double precision (a-h,o-z)
	 dimension al(nrows),  ad(nrows),  ar(nrows)
c	decomposes a tridiagonal matrix a into l*u
c	l  lower triangular (tridiagonal) with all ones on diagonal
c	u  upper triangular (tridiagonal)
c	when calling tridec should have  al(i)=a(i,i-1)  ad(i)=a(i,i)
c	       ar(i)=a(i,i+1)
c	u  is put into ad,  ar
c	l  is put into al  (one knows that the diagonal of  l  is  1.0)
c	al(1)  ar(nrows)   are not used in tridec or in trisol
	 if(nrows  .le.  1)  return
	           nrm1=nrows-1
	           do 40  i=1,nrm1
	 ip1=i+1
c	add  xi*row  i    to row (i+1)   which cancels out a(i+1,i)  and
c	      add  xi*a(i,i+1)  to  a(i+1,i+1),     xi=a(i+1,i)/a(i,i)
	 xi=-al(ip1)/ad(i)
	 ad(ip1)=ad(ip1)+xi*ar(i)
40              al(ip1)=-xi
c	statement  40   gets  l,  l(i+1,i)=-xi
	 return
	 end
c=======================================================
c=======================================================
	 subroutine trisol(al,ad,ar,nrows,rtside,soln)
	 implicit double precision (a-h,o-z)
	 dimension al(nrows),  ad(nrows),  ar(nrows)
	 dimension rtside(nrows),  soln(nrows)
c	tridiagonal backsolve, decomposed matrix from tridec is in al,ad,ar
c	solving l*u*soln=rtside
c	rtside is destroyed in trisol
c	al(1)   ar(nrows)   are not used in tridec or in trisol
c	first solve   l*(z=u*soln)   =   rtside
c	then solve   u*soln  =  z
	 if(nrows  .eq.  1)  soln(1)=rtside(1)/ad(1)
	 if(nrows  .eq.  1)  return
	 nrm1=nrows-1
c	solve  l*z=rtside,  put answer in rtside
c	 rtside(1)=rtside(1)
	 do  20  i=2,nrows
20              rtside(i)= rtside(i)- rtside(i-1)*al(i)
c	now solve  u*soln= (z=rtside)
	 soln(nrows)=rtside(nrows)/ad(nrows)
	 do  40  i=1,nrm1
	 nmi=nrows-i
40            soln(nmi)=(rtside(nmi)-soln(nmi+1)*ar(nmi))/ad(nmi)
	 return
	 end
	 
	 subroutine setup
c	set  up  difference  scheme coefficients
	 parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
 	 implicit double precision (a-h,o-z)
	 common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	 character chartin*1,infile*20,outfile*20
	common/comc/chartin,infile,outfile
	common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping
	 kdt   =  0
	 t  =  zero
	 xrt  =  zero
	 npts  =  0
	 xpta(npts)  =  zero
	 ua(npts)  =  zero
	 x  =  zero
	 if(initu  .eq.  1)  ua(npts)  =  truesol(x,t)
	 nint  =  0
	 ktop  =  0
	 kc =  0
	 do  100  ilayer = 1,  nlayers
	       scalea(ilayer) = solwpa(ilayer) * rhoa(ilayer) / 100.d0
	       n = nintlra(ilayer)
	       ktop = ktop + n + 1
	       xrt = xrt + thicka(ilayer)
	       dx = thicka(ilayer)/dble(n)
c	       left  endpoint
	       kc = kc + 1
	       xa(kc) = x
c	       get finite difference grid in layer ilayer
	       do  50  isublr = 1,  n
c	              right endpoint
	              npts = npts + 1
	              nint = nint + 1
	              x = xpta(npts-1) + dx
	              xpta(npts) = x
	              dxa(nint) = dx
	              ua(npts) = zero
	              rtside(npts) = zero
	              if(initu  .eq.  1)  ua(npts)  =  truesol(x,t)
	              kc = kc + 1
	              xa(kc) = x
50                    continue
100           continue
c	point indices started at 0, to get actual number
c	of points add 1 to npts
	 npts = npts + 1
	 celeft = scalea(1) * outrhl / 80.d0
	 cert = scalea(nlayers) * outrhr / 80.d0
	 nptsm1 = npts - 1
	 nptsm2 = npts - 2
	 nunkn = npts - 2
	 if(isymrt  .eq.  1)  nunkn = npts - 1
	 nunknt = nint - 1 + isymrt
	 if(nunknt  .ne.  nunkn)  stop'nunknt ne nunkn in setup'
	 if(kc  .ne.  ktop)  stop'kc ne stop in setup' 
c	get toutfac which converts input value tf of final time in units
c	specified by chartin into second,  i.e.,  tfinal  (sec) = tf * toutfac
	 if(chartin  .eq.  's')  toutfac = one 
	 if(chartin  .eq.  'i')  toutfac = 60.d0
	 if(chartin  .eq.  'h') toutfac = 60.d0 * 60.d0
	 if(chartin  .eq.  'd') toutfac = 60.d0 * 60.d0 * 24.d0
	 if(chartin  .eq.  'w') toutfac = 60.d0 * 60.d0 * 24.d0 * 7.0d0
	 if(chartin  .eq.  'm') toutfac = 60.d0 * 60.d0 * 24.d0 * 30.43667d0
	 if(chartin  .eq.  'y') toutfac = 60.d0 * 60.d0 * 24.d0 * 365.24d0
c	get  tfinal  (seconds)  and t  (second)
	 tfinal = tf * toutfac
	 dt = tfinal / dble(ntsteps)
	 tout = t / toutfac
c	get  finite difference scheme matrix
c********* at present assume isymrt = 0 *****************
	 if(isymrt.ne.0)stop'  isymrt ne 0 have not yet programmed case'
	 onema = one - alpha
	 twodt = two * dt
	 sixdt = 6.0d0 * dt
	 threedt = 3.0d0*dt
	 i = 0
	 do 300 ilayer = 1, nlayers
	        n = nintlra(ilayer)
	        nn = n
	        if(n  .lt.  2)  stop'nintlra(ilayer)  lt 2 in setup'
	        if(ilayer  .eq.  nlayers) nn = n-1
	        do 200 isubint = 1, nn
	             i = i + 1
c	
c	    x(i-1)        interval  i          x(i)           interval  i+1         x(i+1)
c	         si = scalea                          sip = scalea
c	        ei = scalea*d                     eip = scalea*d
c	
	        si = scalea(ilayer)
	        sip = scalea(ilayer)
	        ei = scalea(ilayer) * diffcfa(ilayer)
	        eip = scalea(ilayer) * diffcfa(ilayer)
	        if(isubint  .eq.  n) sip = scalea(ilayer+1)
	if(isubint  .eq.  n) eip = scalea(ilayer+1) * diffcfa(ilayer+1)
	        dxi = dxa(i)
	        dxip = dxa(i+1)
c	       dirichlet data at left endpoint so unknown number k is at
c	       right endpoint of interval number k
	        if(lumping  .eq.  1)  then
c	mass matrix lumped
	        al(i) = -ei*alpha/dxi
	ad(i) = (si*dxi + sip*dxip)/twodt  +
     &                                  (ei/dxi + eip/dxip)*alpha
	        ar(i) = -eip*alpha/dxip
c	 for  forming rtside
	        aln(i) = ei*onema/dxi
	adn(i) = (si*dxi + sip*dxip)/twodt  -
     &                                   (ei/dxi + eip/dxip)*onema
	        arn(i) = eip*onema/dxip
	 else
c	mass matrix not lumped
	        al(i) = -ei*alpha/dxi + si*dxi/sixdt
	ad(i) = (si*dxi + sip*dxip)/threedt  +
     &                                    (ei/dxi + eip/dxip)*alpha
	ar(i) = -eip*alpha/dxip + sip*dxip/sixdt 
c	 for  forming rtside
	        aln(i) = ei*onema/dxi + si*dxi/sixdt
	        adn(i) = (si*dxi + sip*dxip)/threedt  -
     &                                  (ei/dxi + eip/dxip)*onema
             arn(i) = eip*onema/dxip + sip*dxip/sixdt
	 end if
200                 continue
300            continue
	 if(i  .ne.  nunkn)  stop'i ne nunkn in setup'
	 call tridec(al,ad,ar,nunkn)
	 return
	 end

cc      FORTRAN SOURCE FILE : outp.for
	 subroutine outp(npppsi)
c	do printout of results 
	 parameter(nlmax=20,nptsmax=800,kmax=nptsmax+nlmax)
 	  implicit double precision (a-h,o-z)
	 common/comr/rhsol,xrt,thicka(nlmax),
     & diffcfa(nlmax),rhoa(nlmax),outrhl,outrhr,celeft,cert,
     &scalea(nlmax),solwpa(nlmax),xpta(0:nptsmax),
     &dxa(nptsmax),ua(0:nptsmax),ca(kmax),xa(kmax),wtpera(kmax),
     &h2olra(nlmax),totalw,t,tfinal,tf,tout,toutfac,dt,alpha,
     &al(nptsmax),ad(nptsmax),ar(nptsmax),rtside(nptsmax),
     &aln(nptsmax),adn(nptsmax),arn(nptsmax),scra(nptsmax),
     &zero,one,two
	 character chartin*1,infile*20,outfile*20
	common/comc/chartin,infile,outfile
	common/comi/nlayers,nintlra(nlmax),nint,npts,nptsm1,
     & nptsm2,nunkn,ktop,ntsteps,kdt,levt,initu,krtbc,kleftbc,
     & isymrt,lumping
	 emax = zero
	 write(66,600)
	 write(66,600)
	 write(67,600)
	 write(67,600)
600            format('            ')
	 dtout = dt/toutfac
	 if(chartin  .eq.  's')  write(66,6051)   kdt,  tout,  dtout 
	 if(chartin  .eq.  'i')   write(66,6052)   kdt,  tout,  dtout
	 if(chartin  .eq.  'h')  write(66,6053)   kdt,  tout,  dtout
	 if(chartin  .eq.  'd')  write(66,6054)   kdt,  tout,  dtout
	 if(chartin  .eq.  'w') write(66,6055)   kdt,  tout,  dtout
	 if(chartin  .eq.  'm') write(66,6056)   kdt,  tout,  dtout
	 if(chartin  .eq.  'y')  write(66,6057)   kdt,  tout,  dtout
	 if(chartin  .eq.  's')  write(67,6051)   kdt,  tout,  dtout 
	 if(chartin  .eq.  'i')   write(67,6052)   kdt,  tout,  dtout
	 if(chartin  .eq.  'h')  write(67,6053)   kdt,  tout,  dtout
	 if(chartin  .eq.  'd')  write(67,6054)   kdt,  tout,  dtout
	 if(chartin  .eq.  'w') write(67,6055)   kdt,  tout,  dtout
	 if(chartin  .eq.  'm') write(67,6056)   kdt,  tout,  dtout
	 if(chartin  .eq.  'y')  write(67,6057)   kdt,  tout,  dtout
6051         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  seconds      dt = ' ,e11.3, '  seconds')
6052         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  minutes      dt = ' ,e11.3, '  minutes')
6053         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  hours      dt = ' ,e11.3, '  hours')
6054         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  days      dt = ' ,e11.3, '  days')
6055         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  weeks      dt = ' ,e11.3, '  weeks')
6056         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  months      dt = ' ,e11.3, '  months')
6057         format('  time step # = ' ,i6, '        time = ' ,1pe11.3,
     &               '  years      dt = ' ,e11.3, '  years')
c
	 write(66,6061)
	 write(67,6061)
6061     format(' u = scaled solution = c/scale factor, c =',
     &               'H2O concentration  (g/cc)')
	 write(66,6062)
	 write(67,6062)
6062   format('wtper = (g H2O)/(100 g  dry material) =(100 c)/(rho)')
	 write(66,6063)
	 write(67,6063)
6063         format('  scale factor = solwpa * rho / 100  ')
c	
c	get concentration array ca, and weight percent array wtpera
	 ipt = 0
	 kc  = 0
	 do 99 ilayer = 1, nlayers
	                n = nintlra(ilayer)
c	                left endpoint
	                kc = kc + 1
	                ca(kc) = ua(ipt) * scalea(ilayer)
	                wtpera(kc) = ca(kc) * 100.d0  /  rhoa(ilayer)
	      do 49 isublr = 1, n
c	                     right endpoint
	                     ipt = ipt + 1
	                     kc = kc + 1
	                     ca(kc) = ua(ipt) * scalea(ilayer)
	                     wtpera(kc) = ca(kc) * 100.d0 /rhoa(ilayer)
49                   continue
99              continue
c	now get h2olra, totalw
	 half = 0.5d0
	 kc = 0
	 totalw = zero
	 do 101 ilayer = 1, nlayers
	                 n = nintlra(ilayer)
	delta = thicka(ilayer)/dble(n)
	scrh = zero
c	left endpoint
	kc = kc + 1
	       do 51 isublr = 1, n
c	     left endpoint
	     scrh = scrh + half * ca(kc)
c	     right endpoint
	     kc = kc + 1
	     scrh = scrh + half * ca(kc)
51                   continue
	h2olra(ilayer) = scrh*delta
	totalw = totalw + scrh*delta
101              continue
	  itop = 0
	  do  100 ilayers = 1,nlayers
	                 n = nintlra(ilayer)
	levx = n/npppsi
	if(levx  .lt.  1)  levx = 1
	ibot = itop
	itop = ibot + n
	write(66,610)  (xpta(k),  k=ibot, itop, levx)
610                  format('    x =  ', 1pe14.4,4e14.4)
	write(66,611)  (ua(k),  k=ibot, itop, levx)
611                  format('    u =  ', 1pe14.4,4e14.4)
	if(initu  .eq.  1)  then
	           do  50 k = ibot,  itop
	 xx = xpta(k)
	 scra(k+1) = ua(k) - truesol(xx, t)
	 scre = abs(scra(k+1))
	 if(scre  .gt.  emax)  emax = scre
50                        continue
	    write(66,610)  (scra(k+1),  k=ibot, itop, levx)
612                     format('    e =  ', 1pe14.4,4e14.4)
	                 end if
	write(66,600)
100           continue
	           if(initu  .eq.  1)  write(66,600)
	           if(initu  .eq.  1)  write(66,613)  emax
613           format('  maximum error =  ', 1pe14.4)
	           if(initu  .eq.  1)  write(66,600)
	           if(initu  .eq.  1)  write(66,600) 
c	
	           if(initu  .eq.  0)  then 
	write(66,600)
	write(66,600)
	write(67,600)
	           itop = 0
	 do  110 ilayers = 1,nlayers
	                 n = nintlra(ilayer)
	levx = n/npppsi
	if(levx  .lt.  1)  levx = 1
	ibot = itop + 1
	itop = ibot + n
	write(66,710)  (xa(k),  k=ibot, itop, levx)
710                  format('    x =  ', 1pe14.4,4e14.4)
	write(66,711)  (ca(k),  k=ibot, itop, levx)
711                  format('    c =  ', 1pe14.4,4e14.4)
	write(66,712)  (wtpera(k),  k=ibot, itop, levx)
712                  format('    wtper  =  ', 1pe14.4,4e14.4)
	write(66,600)
110          continue
	itop0 = 0
	itop1 = 0
	write(67,6767)
6767  format('  x =  location (cm)   u = called conc.  H2O   conc.',
     &              '  (g/cc)       g  H2O/100g  material')
	write(67,600)
	 do  120 ilayers = 1,nlayers
	                 n = nintlra(ilayer)
	levx = n/npppsi
	if(levx  .lt.  1)  levx = 1
	ibot0 = itop0
	itop0 = ibot0
	ibot1 = itop1 + 1
	itop1 = ibot1 + n
	k1 = ibot1
	do 119  k=ibot0,itop0,levx
	write(67,6768)  xpta(k),ua(k),ca(k1),wtpera(k1)
6768                format(1pe14.4,e20.4,e19.4,e23.4)
	k1 = k1 + levx
119              continue
120              continue
	          end if
	write(66,600)
	write(66,600)
	write(67,600)
	write(67,600)
	write(66,749)
	write(66,749)
749    format('  water in each layer as gram in a 1 cm*cm cross',
     &      '  section of the material ')
	write(66,750)  (k,  h2olra(k),  k = 1,nlayers)
	write(67,750)  (k,  h2olra(k),  k = 1,nlayers)
750   format(2(' layer =',i3,' H2O  in layer =  ', 1pe11.3,'  '))
	write(66,600)
	write(66,600)
	write(67,600)
	write(66,751)
	write(67,751)
751   format('  total amount of water (gram) in a 1 cm*cm cross',
     &  '  section of the full material ')
	write(66,752) totalw
	write(67,752) totalw
752         format('  total amount of H2O =  ', 1pe14.4)
	write(66,600)
	write(67,600)
	          return
	          end









c GLOSSARY OF VARIABLE NAME IN THE PROGRAM diffuse
c EXPLANATION OF THE ALGORITHM USED IN THE PROGRAM diffuse
c	THE PROGRAM DIFFUSE calculates the time history of the water
c	concentration in a sample composed of distinct layers having different 
c	material properties (which are constant inside each layer). The
c	sample I considered to be equivalent to a cylinder exposed to a
c	constant relative humidity at it ends, and insulated along its
c	lateral surface, with its axis along the x direction starting at x=0.
c	
c	As x varies from 0 to x=xrt (the right endpoint of the cylinder), one
c	passes through the various layers in the sample. For a given fixed
c	value of x, there is no change in the material properties or in the
c	H2O concentration a y and z vary, so the mathematical model only
c	involves 1 space dimension
c	
c	ci   designates that the variable being described is an output 
c	          variable (read in from a data file at the start of the program)
c*****     NOTE, except for comment which I rea in with a72 format,
c*****     ALL DATA REA IN FROM THE INPUT DATA FILE, infile, IS
c*****     READ IN USING FREE FORMAT.
c	co   designates that the variable being described is an output
c	           variable (whose value is computed by DIFFUSE).
c	cp   designates that the variable being described is a parameter
c	           constant in the program (e.g.,  a maximum for the number of
c	           finite difference points) which can be changed (in the include
c	           file comincl) if necessary.
c	cs   designates that the variable being described is an input
c	          variable which is read from the terminal at the start
c	          of the program (the names of the input an output files).
c	variables without an i  o  p  or   s  designator are internal variable used 
c	by DIFFUSE.
c	
c	DIFFUSE used an include file,  comincl, containing various
c	                  parameter, common block and declaration statements
c	                  neee by the routines in DIFFUSE.
c	
c	CHARACTER VARIABLES
c i	chartin is a character variable of length 1, specifying time units,
c i	comment is a character variable of length 72, used to input and
c	                output comments on a given run,
c	cscr is a character scratch variable of length 3,
cs	infile is a character variable of length 20 used for the input data
c	         file name (infile is read in from the terminal at the start
c	         of the program-follow instructions printed on screen),
c	outcol  is a character variable of length 23 used for the file name
c	            of the output file containing the output data in column format,
c	            outcol = <outfile with trailing blanks delete>col,
c****  when typing in infile and outfile on the keyboard, do NOT use
c****  any blank spaces, do NOT enclose the names in apostrophes or
c****  quotes, and o not have the file name exceed the
c****  maximum file name length on the computer being used.
c****  the name of the file containing output in column format is
c****  defined to be <the output file name type onto the screen>col
c****  which thus has 3 more characters in its name than does the 
c****  regular output file.
c	
c	all other variables starting with the letter a-h,o-z are
c	double precision,
c	all other variables starting with the letter i-n are integers.
c	
c******      UNITS ABBREVIATIONS
c	
c	g = gram
c	sec = second
c	cm = centimeter
c	cc = cubic cm
c	
c******      PHYSICAL VARIABLES
c	
c	xrt:  the sample (axis of the cylinder) is considered to be located
c	       between  x=0  and x=xrt  (cm).
ci	nlayers:  number of layers of (different) materials in the ample
c	              under consideration.
cp	nlmax: maximum number of layer allowed  (this is a programming
c	            parameter which can be increased if necessary; its current
c	            value in the include file comincl   is 20).
ci	thicka(ilayer):  array (vector) of thickness (in cm)
c	                        of each layer (ilayer = 1,,nlayers) in the sample.
ci	diffcfa(ilayer):  array (vector) of diffusion coefficients in
c	                         each layer (ilayer = 1,,nlayers) (units are
c	                         cm*cm/sec).
ci	rhoa(ilayer): is the array (vector) of density values of ry material
c	                        for each layer (ilayer = 1,,nlayers)
c	                        (units are g/cc).
c	isymrt: if 1 then the right endpoint I seale (no moisture diffusion
c	takes place through the right endpoint of the sample). The
c	boundary condition at the right endpoint is then the no flux
c	(symmetry) boundary condition.
c	If isymrt I not 1, then the right endpoint is
c	exposed to the ambient moisture. The latter is ALWAYS 
c	the case for the current version of DIFFUSE.
ci	outrhl: the relative humidity (in percent) outside of the 
c	            left endpoint (x=0) of the sample.
ci	outrhr: the relative humidity (in percent) outside of the 
c	            right endpoint (x=xrt) of the sample.
c	cross section: for purpose of  1-dimensional diffusion calculations
c	we consider the situation of a 1 square cm cross                       
c	section of material, so a length of 1 cm corresponds to
c	a volume of 1 cc (cubic cm). The numerical solution
c	procedure is a mass conservative finite difference
c	method. When alpha = 1. and lumping = 1 (see below):
c	this method satisfies a maximum
c	principle (the range of the solution values
c	must stay within that of the initial and boundary
c	data), and it is equivalent to the
c	lumped fully implicit finite element method
c	with linear elements for this diffusion problem.
c	concentration: denoted by  c  or   c(x)   or  c(x,t), is the concentration 
c	of H2O at a particular location  x  in the sample,
c	at some time  t  (in units of grams of  H2O/cc,  or
c	grams of H2O/100 grams of dry materials).
c	The solution of the diffusion equation is done  in
c	terms of grams of H2O/cc, and UNLESS OTHERWISE NOTED,
c	ALL CONCENTRATIONS ARE IN GRAMS OF H2O/cc.
c	celeft:  equilibrium concentration of H2O in the material at the left
c	            endpoint of the sample when maintained at the relative
c	            humidity value outrhl.
c	cert:  equilibrium concentration of H2O in the material at the right
c	            endpoint of the sample when maintained at the relative
c	            humidity value outrhr.
c	scalea(ilayer): concentration scale factor array (vector)
c	(ilayer = 1,,nlayer).  If the material in layer
c	ilayer and some "standard" material are both
c	maintained to equilibrium at some given relative
c	humidity, then the H2O concentration in the material
c	in layer ilayer will equal calea(ilayer)  * the H2O
c	concentration in the "standard" material. Thus at
c	equilibrium the normalized or scaled concentrations
c	           c(ilayer)/scalea(ilayer)
c	will all be the same.  It is being ASSUMED that right
c	at the interface   I  between 2  materials
c	
c                             material  i                   material i+1
c                        ---------------------- -------------------------
c                                                     I 
c	
c	
c	
c	
c	the relationship   c(i)/scalea(i) = c(i+1)/scalea(i+1)
c	always holds  (even before equilibrium is attained).
c	scalea(ilayer)  is defined as a dimensionless constant.
ci	rhsol:   relative humidity (%) t which equilibrium moisture
c	           solubility data (solwpa) for all the materials have been
c	           measured.
ci	solwpa(ilayer): is the array (vector) of the equilibrium solubility
c	of H2O I the material of layer ilayer
c	(ilayer= 1,,nlayers) maintained at rhsol% relative
c	humidity. units are g H2O / 100 g dry material =
c	weight percent of H2O.
c***********       NOTE immediately after solwpa has been read in and
c	thn printed out,  it is multiplied by
c	                 80.  /  rhsol
c	to rescale it to correspond to saturation levels
c	at 80% relative humidity (since that is what
c	was assumed for the original version of this code).
c	NOTE also celeft = (calea(1) * outrhl  /  80.)*1  g/cc
c	and   cert  =  (scalea(nlayers) * outrhr  /  80.)*1  g/cc
c	
c****** NUMERICAL SOLUTION PROCEDURE VARIABLES / OUTPUT VARIABLES
c	
ci	nintlra(ilayer): array(vector) of number of subinterval each layer
c	(ilayer = 1,,nlayers)  is subdivided into for the
c	finite difference method.  Thus the thickness of each
c	subinterval in layer ilayer is
c	thicka(ilayer) / nintlra(ilayer).
c****                     the value of nintlra(ilayer) MUST be at least 2
c****                     for each material layer.
c	nint:  total number of subintervals into which  (0,  xrt) is
c	         subdivided for the numerical solution procedure. nint is the sum
c	         over ilayer = 1,,nlayers of nintlra(ilayer).
c	npts:  is the number of computational points where  c(x)  is calculated,
c	          including the two boundary points where the value of  c(x)  is
c	          specified;  npts = nint + 1
cp	nptsmax: maximum number of computational points allowed (this is a
c	                programming parameter which can be increased if necessary;
c	                its current value in the include file comincl   is 800).
c	nptsm1: is npts - 1
c	nptsm2: is npts - 2
c	nunkn:  is the number of computational points where the unknown value
c	             of  c(x)  is calculated at each time step = npts - 2 = all the
c	             grid point - the two ebdpoints where (known) boundary data
c	             for  c  is specified.  (When isymrt = 1, nunkn = npts -  1)
c	xpta(ipt): array (vector) of  x  point values at which  c(x)  is
c	                calculated  (ipt = 0,,nptsm1)  (units for xpta is cm).
c	                The two endpoints  (x=0 and x=xrt)  are inclue in xpta.
c	dxa (int):  is the mesh spacing array  (interval size array);
c	                 dxa(int) = xpta(int)-xpta(int-1)   (for  int = 1,,nint).
co	   ua(ipt):  is the array (vector) of the approximate values at a given
c	                time  t  for   c(xpta(ipt),t)/scalea(ilayer) where ilayer is the
c	number of the material layer containing   xpta(ipt)   
c	(ipt = 0,,npts-1).  Note from the discussion of scalea,
c	there I no ambiguity in the value of ua at the interface
c	between two layer of material (which is the rationale for
c	the definition of  u = c/scalea).  THE NUMERICAL DISCRETIZATION
c	(APPROXIMATION) OF THE DIFFUSION PROBLEM IS FORMULATED
c	AND SOLVED IN TERM OF ua.
c	Note the boundary conditions for   u  are  (cf.  the specific
c	definition of scalea)  u  =  outrhl / 80.  at  x=0,
c	an  u = outrhr / 80. at  x=xrt.
co	ca(k),  xa,  ktop:  ca(k) is an array (vector) of H2O concentration
c	values in the sample  (k = 1,,ktop).  The  x  value
c	corresponding to ca(k)  is  xa(k).  The values of xa 
c	range over the values of xpta. At a point
c	xa(k) = wpta(ipt) which is not an interface between
c	two layers,
c	ca(k) = ua(ipt) * scalea(ilayer)
c	where ilayer is the index of the layer containing
c	xpta(ipt).  AT INTERFACE POINTS BETWEEN TWO MATERIAL
c	LAYERS THERE ARE 2 VALUES OF  ca  GIVEN;
c	                   ca(k) = ua * (scala value from left side of interface)
c	               ca(k+1) = ua * (scala value from right side of interface)
c	               in which case  xa(k) = xa(k+1) = the interface location.
cp	kmax:  maximum value of ktop allowed  (this is a programming
c	parameter which can be increased if necessary-- it is currently
c	set to nptsmax+nlmax in comincl which will always be adequately
c	large).
co	wtpera(k):  is the array (vector) of concentration values in units of
c	g  H2O / 100 g dry material   (weight percent of H2O)
c	for  k = 1,, ktop.  wtpera(k) =
c	(ca(k) / rhoa(ilayer corresponding to  xa(k)) * 100.
c	At a pair of indices k, k+1 corresponding to an interface,
c	rhoa in the denominator takes the value from the left,
c	then the right side of the interface at  k, k+1
c	respectively.
c	e.g.,  a concentration weight percent of  .6 means  .6  g  H2O/100  g dry
c	material.
c	h2olra(ilayer): array (vector) of total amount of H2O  (grams) in
c	each layer  (ilayer = 1,,nlayers).  This is 
c	obtained from  ca  and  xa (trapezoidal quadrature).
c	The cross section size is 1 square cm.
co	totalw:   is the total amount of water in the sample (g  H2O) = the sum
c	over ilayer = 1,,nlayers of h2olra(ilayer).
c	The cross section size is 1 square cm.
c	t:  is the current time value at which the approximate solution of the
c	    diffusion problem is being obtained (unit are seconds).
c	tfinal:  final time value for which the approximate solution is to be
c	           obtained (unit are seconds).
ci	chartin:  character variable length 1 which specifies the unit for
c	the input value of the final time,  tf,  and the units 
c	for values of time,  tout,  in the output file.
c	chartin = 's' for seconds,  'i' for minutes, 'h' for hours,
c	'd' for days,  'w' for weeks,
c	'm' for months (=30.43667 days),
c	or  'y' for years (=365.24)
ci	tf:  input value of the final time at which the concentration in the 
c	     sample is to be obtained.  (Units for  tf  and tout are specified
c	     by chartin.)
co	tout:  variable containing current ouput time value.
c	toutfac:  tout = t / toutfac  (convert time in sec to output time
c	              units specified by chartin).
ci	ntsteps:  total number of time step to be use to obtain the
c	              approximate solution of the diffusion equation up to
c	              t = tfinal  (time steps are kdt = 1,2,,nsteps).
c	kdt: time step number  (kdt = 1,2,,nsteps) (kdt = 0 during
c	        input of data and initialization of the numerical solution 
c	        procedure).
c	dt:  time step size  (seconds) = tfinal / ntsteps
ci	levt:  print out approximate solution values every  levt  time steps,
c	         i.e., at  t=0 (initial data),  t = levt*dt, t = 2 * levt*dt, etc.
ci	initu :  initialization indicator, set = 0 to set ua=0  at  t=0;  set = 1
c	           to use the function truesol to initialize  ua (for use in
c	           testing the program with a known analytical solution function)
ci	alpha :  time averaging parameter ; alpha = 1. for pure implicit,
c	             alpha = 5. for Crank-Nicolson (Crank-Nicolson is second order
c	             accurate but theoretically not as stable as pure implicit).
c***            for  this application it is probably best to choose alpha=5.
ci	lumping:  if lumping = 1 then "lump the mass matrix", i.e., do not
c	use the to grid points adjacent to the center point to
c	average the time differencing -- when  alpha = 1.,  lumping=1
c	enforces the discrete maximum principle.  When lumping is 
c	not equal  1 the resulting difference scheme is the 
c	ful finite element method with linear elements.
c***            for  this application it is probably best to choose lumping=0
c	
c	
c*********    COMMENT ON CHOOSING SPATIAL MESH AN TIME STEP SIZE
c	
c	In order to verify that sufficiently many finite difference
c	subintervals (nintlra) have been ue for each material
c	layer, and that a sufficiently small
c	time step size has been chosen (ntsteps large enough)  in order to
c	obtain an accurate approximate solution to the moisture diffusion 
c	problem using DIFFUSE, it is suggested that one do the following.
c	(1) do a printout of the numerical solution values at an "early"
c	time  (well before the first time value one is specifically
c	interested in), and make sure there are no unnatural oscillations
c	in the solution values -- if there are then increase ntstep.
c	(2) double ntsteps and double each nintlra(ilayer) value and check
c	to see if the changes in the numerical solution values at the times
c	and spatial locations in which one is interested have changed
c	by more than an acceptable amount. Keep increasing the resolution
c	of the spatial grid and time stepping until the numerical results
c	vary by les than whatever tolerance is desired.
c	
c	Note that one must have a separate "numerical scheme layer" for each
c	distinct material layer. However one is free to have adjacent
c	specified material layers whose physical properties are identical.
c	For example in a sample with same material from  x=2  to  x=4,
c	one is free to use one layer from x=2  to  x=4, or one could use 
c	two layers (x=2  to 3, and  x=3  to  4) or
c	one could take  x=2  to  2.5  an x=2.5  to  4  etc. One can use
c	this to, e.g., "grade the numerical mesh" within one material, or vary
c	the amount of printout in one layer without changing the Fortran.
c	
ci	comment:   up to five lines of comments  (72 characters each) at the
c	end of the input file (after the line where the value of
c	the variable lumping is specified) are read into the
c	character variable comment (with a72 format)
c	and are printed out at the beginning of the output file.
c	npppsi:   governs the number of point printed (in the output file)
c	for each material layer. At present npppsi is 4 which means
c	that, a long a npppsi evenly divides nintlra(ilayer),
c	there will be 5 output points from material layer
c	number ilayer (the 2 endpoints and 3 evenly spaced interior
c	points). the value of nppsi is set at the beginning of 
c	the main program DIFFUSE.
c	kleftbc,  krtbc: indicate type of boundary data at the left, right
c	                        endpoints -- not used within the current version
c	                        of DIFFUSE.
c	al, ad, ar, rtside: vectors used in the finite difference method 
c	                           solution procedure for evaluating  ua  at each
c	                           time step.
c	aln, and, arn:  vectors used in forming rtside.
c	scra:  is a scratch array.
c	zero: the constant 0.
c	one: the constant 1.
c	two: the constant 2.
c	
c	
c******    FORMULATION OF THE DIFFUSION PROBLEM IN TERMS OF   u=c/scalea
c	
c	The flux  F at a given location is -d(x)c'(x)  where  c' denotes the
c	partial derivative of  c with respect to  x. At an interface point
c	the limits of this expression as  x  approaches the interface from
c	either side must be the same.  Setting  e(x) = d(x)*scalea(ilayer)
c	where ilayer is the layer containing the point  x, the flux F
c	is then  -d(x)*scalea * c'(x)/scalea = -e(x)u' (x), and as decribed
c	above,  u(x) is continuous across the material interfaces.
c	We now consider a mass balance finite difference equation at
c	the finite difference grid point  x(i) in the following diagram:
c	
c	
c            u(i-1)            F(i)          u(i)             F(i+1)            u(i+1)
c               o                                  o                                          o
c            x(i-1)        x(i-1/2)        x(i)            x(i+1/2)          x(i+1)
c                                 e(i)                                e(i+1)
c                             scalea(i)                         scalea(i+1)
c	
c	
c	
c	
c	
c	
c	Note in the diagram we are using the notation scalea(i) = the value 
c	of scalea(ilayer) for the material layer containing  x(i-1/2).
c	
c	The net amount of H2O passing into the interval [x(i-1/2), x(i+1/2)]
c	during an interval of time dt  is    dt*(-F(i+1)+F(i)) = (approximately)
c	              dt * e(i+1) * [(u(i+1) - u(i)) / (x(i+1) - x(i))]
c	           -  dt * e(i) * [(u(i) - u(i-1)) / (x(i) - x(i-1))]
c	(To make the unit come out correctly, one would multiply by the
c	1 square cm cross section of the cylinder.)
c	by conservation of mass, this is balanced by the change deltau (during
c	the time dt) in the scaled concentration  u  in the interval
c	[x(i-1/2), x(i+1/2)]  (i.e., in the volume formed by sweeping out the
c	1 square cm cross section over this interval). The change deltau
c	result in the (approximate) change of H2O content in the interval
c	given by
c	                   .5*(x(i+1) - x(i)) * scalea(i+1) * deltau
c	                + .5*(x(i) - x(i-1)) * scalea(i) * deltau
c	
c	Setting the net amount of H2O passing into the interval (13-14 lines
c	above) equal to the expression immediately above gives the finite
c	difference equation at the (interior) grid point  x(i). Values for
c	u  at the two endpoints  (x=0,  xrt) are specified (determined by the
c	ambient relative humidity). THESE EQUATIONS ARE IDENTICAL TO THE
c	RESULTS OF APPLYING THE LUMPED FINITE ELEMENT METHOD WITH 
c	LINEAR ELEMENTS TO THE EQUATION :
c	
c	                 scalea (du/dt) = (e  u')'
c	or
c	                 scalea (du/dt) = (scalea*d  u')'
c	or
c	                 (solwpa*rhoa/100.) (du/dt) = ((solwpa*rhoa/100.)*d  u')'
c	
c	where  du/dt designates the partial derivative of  u  with respect to  t.
c	Note the boundary conditions for  u  at the two endpoints are given by
c	                 u = outrhl / 80.  at x = 0, and
c	                 u = outrhr / 80.  at x = xrt.
c	The concentration of H2O is then given by   c = scalea * u.
c	
c	The finite difference equation above result in a tridiagonal 
c	linear system of equations which needs to be solved whenever
c	u  (and thus  c) are obtained at a time  t+dt using the current (known)
c	values at time  t  and the boundary data.
c	
c******    SUBROUTINE IN diffuse
c	
c	diffuse: main program
c	initl: read in data from input file physical constants and
c	        specifications for the finite differences method).
c	setup:  sets up difference scheme coefficients.
c	advance:  advances the approximate solution one  dt  time step.
c	outp:  writes specified output to output file.
c	tridec:  routine used to decompose (factor)  the
c	tridiagonal linear system (matrix) from the
c	finite difference equations for advancing the approximate
c	solution of the diffusion problem one time step.
c	trisol:  routine used at each time step to solve (backsolve) the
c	factored tridiagonal linear systems of finite difference
c	equations in order to get the approximate solution
c	values at the next time level.
c	truesol, source :  functions used to test the program
c	                            using a known exact solution.
c	dirl, dirr : finctions used to get the boundary values (at  x=0,  x=xrt)
c	for the scaled solution  u.
c	The Fortran functions dirl and dirr are functions of
c	time.  In the present version of DIFFUSE they only depend
c	on the constants outrhl, outrhr (their argument,  i.e.,
c	time, is not used). The structure of DIFFUSE is such that
c	time dependent boundary data could be programmed into
c	dirl  and dirr.  Note some compilers may issue a message
c	that the argument (time) of dirl and dirr is not used.
c	In this case, such a message is of no consequence.
c	parameter specifications, variable declarations and common blocks
c****  the file comincl houl reside in the same directory which
c****  contain the Fortran program routines of DIFFUSE.
c***  SAMPLE INPUT DATA (schematic sample)
c***  this is not "real" physical data, but serves to explain the
c***  layout of the input data. See actual sample input data and
c***  corresponding output from running the program on a
c***  real application. 
c***  SAMPLE INPUT DATA (free format expect for a72 format for comment)
c***  sample input data starts on line below.
c 3 nlayers  # material layer; for each layer, read in its physical props.
c 40.  rhsol   rel  humidity at which solwpa data has been measured  (%)
c 1.         thicka:  thickness of this layer in centimeters (cm)
c 1.         diffcfa:  diffusion coeff. of material in this layer (cm*cm/sec)
c 1.         rhoa:  density of the (dry) material in this layer (g/cc)
c 1.         solwpa:  equil.  sol.  of H2O in material in this layer (wt %) @rhsol
c  2.       thicka:  thickness of this layer in cm
c  2.       diffcfa:  diffusion coeff. of material in this layer (cm*cm/sec)
c  2.       rhoa:  density of the (dry) material in this layer (g/cc)
c  2.       solwpa:  equil.  sol.  of H2O in material in this layer (wt %) @rhsol
c    3.     thicka:  thickness of this layer in cm
c    3.     diffcfa:  diffusion coeff. of material in this layer (cm*cm/sec)
c    3.     rhoa:  density of the (dry) material in this layer (g/cc)
c    3.     solwpa:  equil.  sol.  of H2O in material in this layer (wt %) @rhsol
c45.  55.  ourhl,  outrhr: rel. hum.  (%) outside of the left,  rt endpt
c10   nintlra: no.  of finite diff. subintervals  1#  for each layer(layer1)
c20   nintlra: no.  of finite diff. subintervals  1#  for each layer(layer2)
c30   nintlra: no.  of finite diff. subintervals  1#  for each layer(layer3)
c  'h'   chartin: character variable of length 1 giving input time units
c8.   tf: final time (units specified by chartin) calc. carried out to
c400   ntsteps: number of finite difference method time steps to get to tf
c20    levt: print out approx. solution values every levt time steps
c0   initu: if  0, start with 0 concentration, if 1 test with known c(x,t)
c0.5 alpha:  time ave. parameter (.5 is Crank-Nicolson=recommended choice)
c0 lumping: 0 means use the full finite element method=recommended choice
c comment1: any remaining lines (up to a maximum of 5) in the input
c comment2: file are considered to be comments which are printed at the
c comment3: beginning of the output file
c comment4: the comment lines (if any) are read in with  a72 format so
c comment5: no quote marks are to be used to delimit the comment text


